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THE mviSCID PRESSUBE FIELD OH THE TIP OF A SEMI-IHFnjITE WHtG 

AHD ITS APPLIGAriOK TO THE FOEMATION OF A TIP VORTEX 

By G. F. Hall, S. J. Shamroth, H. McDonald and W. E. Briley 
United Technologies Research Center 

SUMMARY 


A method has been developed for determining the aerodynamic loads on the 
tip of an infinitely thin, swept, cambered semi-infinite wing at an angle of 
attack which is operating subsonically in an inviscid medium and is subjected 
to a sinusoidal gust. Under the assumption of linearized aerodynamics, the 
loads on the tip are obtained by si;g)e 2 rposition of the steady aerodynamic 
results for euigle of attack and camber, and the unsteady results for the res- 
ponse to the sinusoidal gust. The assumption that the major effects of the 
tip are confined to a relatively small finite spanwise length permits the 
loading over the inboard portions of the wing to be approximated by known 
results for the infinite wing. The tip region is treated by a numerical 
doublet-lattice method. The leading edge singularity in the loading is 
removed by applying a correction factor due to Lighthill. The near field 
disturbance pressures in the fluid surrounding the tip are obtained by assum- 
ing a dipole representation for the loading on the tip and calculating the 
pressures accordingly. In addition to providing insight into the inviscid flow 
behavior in the vicinity of the tip, the near field pressures are also used to 
drive a reduced form of the Havier-Stokes equations which yield the tip vortex 
formation. 

Although the derivation is valid for the general conditions described 
above and several inviscid cases have been calculated, the combined viscid- 
inviscid analysis has only been applied to determining the pressures and 
examining the vortex rollup in the vicinity of an unswept, uncambered wing 
moving steadily at a Mach number of 0.2 at an angle of attack of 0.1 rad. 
Generally, the pressure behavior is as expected. As the wing surface is ap- 
proached in the normal direction, the pressure magnitude increases, culminat- 
ing at one-half the dipole magnitude at the surface . External to the region 
containing either the wing or the wake, the pressures go continuously and 
monotonically to zero as the wing plane is approached. In the wake, the 
magnitude increases as the wake is approached until, near the wake, an abrupt 
reversal in trend occurs and the magnitude goes to zero. Finally, in the 
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vicinity of the trailing edge, the pressure magnitude increases until, near 
the wing surface, a reversal in trend is noted -with the magnitude decreasing 
and terminating at one-half the dipole strength. This behavior reflects the 
change in boundary condition at the trailing edge. 

The viscous tip flow calculation has shown features expected in the tip 
flow such as the gualitatively proper development of boundary layers on both 
the upper and lower airfoil siirfaces. In addition, application of the vis- 
cous solution leads to the generation of a 'circular* 1^6 flow pattern 
above the airfoil suction- surface . 


2 



IMTRODUCTIOIT 


In recent years, there has been a growing interest in the vortex formed 
due to the pressure eciualization condition at the tips of lifting surfaces. 

This tip vortex, whose presence can often be neglected in isolated -airfoil 
aerodynamics, requires a precise definition of its physical characteristics 
in problems involving interference aerodynamics. A primary example of this 
type of problem is the- intersection of a helicopter rotor blade tip -vortex ■ - 
with the next following blade. This intersection represents a potentially 
serious acoustic problem, (Ref. l) and has been shown, to have some effect on 
rotor performance (Ref. 2). In addition to the direct relevance to cxirrent 
helicopter problems, a more complete knowledge of tip vortex dynamics is 
needed to assist in the design of novel helicopter blade tips for optimum 
performance, and to provide basic information on the behavior of trailing tip 
vortices behind large aircraft, to name only two of many applications. 

In the past, the tip vortex roll-up rate and position have been deter- 
mined either by direct experimentation or by analytical techniques utiliz ing 
Trefftz plane calculations (e.g.. Ref. 3). Experimental procedures are reason- 
able for many steady flow problems, but they are not yet sufficien ly sophis- 
ticated to show all the details of the vortex formation. Trefftz plane analy- 
ses generally predict the position of the tip vortex adequately, but yield 
roli-up rates much too slow. Hall (Ref. 4) has modeled the vortex wake behind 
a fixed wing starting impulsively from rest which predicts both a rate and 
position of the tip vortex more compatible with those observed, but ’this tech- 
nique requires large amounts of conpiuter time and storage. Another approach 
to determining the tip vortex behind a lifting surface is a direct numerical 
assault on the Navier-Stokes equations, but this 'approach also requires large 
amounts of computer time and storage. It should be noted that classical 
lifting surface analyses are incapable of treating the roll-up problem; the 
linearization inherent in these methods is not valid in the tip region. The 
analysis of the tip vortex formation is a highly nonlinear real fluid problem. 
Pigure 1 is a sketch depicting this formation. 

The present analysis is an attempt to provide a simplified model of the 
tip vortex formation problem. A boundary layer approach is used in which the 
vorticity laden fluid contained in the tip vortex originates in the region near 
enough the wing surface to be affected by viscosity. The boundary layer condi- 
tion that the pressure field is unaffected by viscous displacement effects is 
assvuned so that this pressure field can be determined from inviscid potential 
theory in which shed vorticity is modeled by discrete vortex elements . It is 
further assumed that the phenomenon is sufficiently confined to the tip region 
that a semi-infinite wing can be used to determine the pressure field. Since 
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the presence of the vortex requires a loading to exist on the wing, angle of 
attack, camber and unsteadiness must be considered, but thickness can be 
neglected, at least as a first approximation. This models the wing as an 
infinitely thin, curved plate in a subsonic compressible flow containing a 
sinusoidal gust. Once the loading on the wing Is determined, the pressure 
field in the fluid around the wing is determined from a pressure dipole solu- 
tion. The viscous roll-up is then obtained from a reduced form of the Wavier - 
Stokes equation. A sketch of the expected tip vortex formation on a flat plate 
tip is shown in Fig. 2.- A con^arison with Fig. 1 shows that the only basic 
difference in the -hypothesized vortex formation process is that on the wing of 
finite thickness the vortex initiation point may lie aft of the leading edge 
whereas on the infinitely thin wing the initiation point is forced to occur at - 
the tip leading edge. 

The approach to be used in this analysis is broadly based on the infinite 
wing analysis of Adamczyfc (Refs.' 5,6) and the doublet -lattice method of Giesing, 
et al (Eef. 7). In brief, Adamczyk examined the problem of an arbitrary g^lst 
encounter of an infinite swept wing in a subsonic flow. The solution, which 
basically treated the sinusoidal gust response, is used in a Fourier super- 
position method for any arbitrary disturbance, and this technique is employed 
herein' for the unsteady vortex formation. Although Adamczyk 's method is posed 
for an infinite wing,, it is known that with appropriate modifications in the 
vicinity of a finite tip, such a two-dimensional method will still be valid at 
some distance inboard of the tip. Hence, the Giesing method is applied over 
the actual tip region and is matched to the Adamczyk solution inboard of the 
tip. The match point location -was determined by numerical experimentation 
and was chosen to yield“a‘ smooth transition between inboard and outboard solu- 
tions. This procedTire is described in the text below. 
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LIST OF SYMBOLS 


A 

b 


c 


C 


P 



D 

E(K) 

F,f 


hi, h2, 


k 




Speed of sound, m/sec 
Influence coefficient 
Semichord, normal to span, in 
Chord, parallel to free stream, m 
Pressure coefficient (p-Po) 

Pressure coefficient change across lifting surface, 

(p(x2-)_p(x2+))2/ppU2 

ACp due to an infinite wing in a given flow 

Sublayer damping factor 

Defined by Eg. 11 

Conplex Fresnel Integral 

Equation of the lifting surface 

Zeroth and first order Bessel functions, respectively 
Wing thickness, m 
metric coefficients 

Zeroth and first order Hankel functions of the first kind 
Gust wave number, 

Gust wave number, x-conponen.t, 

Gust wave number, y-conponent, 

Gust wave number in streamwise direction (Ref. 7), 211 "^ 
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LIST OF SB©OLS (Cont’d) 


k 

L 

P 

Po 

R 

t 

T 

T 

L 

u,v,w 

U 

U 


^2 


w 

wg 

x,y 


x,y,z 


'Xl, Xg 


Xg 
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k/b 

Spanwise length outboard from which the tip effects are 
assumed concentrated, m 

Mixiaag length, m 

Pressure, n/m^ 

Inviscid pressure, n/m^ 

Free stream static pressure, n/m^ 

Radius, m or gas constant 

Time , se c 

Teirperature , ° K 

Lift transfer function 

Velocity conponents in viscous analysis, m/sec 

Freestream velocity, m/sec 

Velocity vector, m/see 

CoHplex gust aijplitude (Ref. 5)j m/see . 

CoBplex gust amplitude, m/sec 
Gust velocity, m/sec 

System coordinates, mj x streaiuwise^ y normal to x and 
directed inboard. Origin at tip leading- edge . 

Cartesian coordinates in viscous analysis, m 

System coordinate , m; x normal to span, x^ parallel to span, 
directed inboard. Origin at midchord of tip. 

Coordinate normal to plane of wing, m 



LIST OP SYMBOLS (Cont’d) 


a 

3 

Y 

e 


Tl 

0 

P 


Po 

P2 

tp 

C) 

(_) 


Subscripts 

's. 

k, 1 
i,0,n,m 


Gust inclination angle, degree, or difference operator 
Prequency parameter, defined by Eq. (22a) 

Cutoff parameter, defined by Eq. (22b) 

Included angle between a point at which induced velocity is 
calculated and the part at which the inducing load is located, 
deg 

Transformed transverse variable 
Wing sweep angle, deg 
Viscosity, kg/m-sec 
Density, kg/m^ 

Freestream density, kg/m^ 

Wing leading edge, radius, m 
Complex velocity potential, m^/sec 

Complex amplitude of velocity potential, m^/sec, or arbitrary 
dependent variable 

Complex aniplitude 

Vector 

Mean value 

Fluctuating value 

■Identifies points at which boundary condition is satisfied * 
Identifies points at which inducing loads are located 
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THEORETICAL FORMULATION 


It is assumed for the purpose of analysis that the semi -infinite wing is 
at rest in a subsonic compressible stream having an imbedded sinusoidal gust 
which is swept past the wing at the stream speed, U. The pertinent geometry 
of the semi-infinite wing is given in Fig. 3. Of particular interest is the 
spanwise distance, L. This length represents the distance within which the 
tip effect is assumed to he dominant; inboard of this distance the flow and 
the loading is essentially two-dimensional. Both the tip and the normal to 
the gust wave fronts are parallel to the freestream due to restraints imposed 
by the particular version of the doublet-lattice method selected. The wing 
is assumed infinitely thin so that thickness can be neglected, hut the effects 
of camher and angle of attack are included in the analysis, it is fm-ther 
assumed that all disturbances are small enough that linearization is permitted. 
The disturbance velocity potential is thus defined by the standard linearized 
acoustic eq[uation for subsonic compressible flow. 




dH 

ax2 


Ob axat 


I 

OqZ 3t2 


= 0 


( 1 ) 


A unique solution to the equation is obtained by applying the necessary boTind- 
ary conditions (Ref. 5). The linearized pressure field is then obtained from 


_P 

Po 



+ U 




-1 

xJ 


( 2 ) 


The necessary boundary conditions can he stated as follows: 


o 

II 

0. 

x<-ytan0, x>c-ytan0, ally, X 2-0 

(3a) 

p=0 

x=c-ytan0, y>0, X2=0 (Kutta Condition) 

(3h) 

Dt 

x>-yton0, x<c-yton0, X 2 = 0, y>0 

(3c) 

(Flow Tangency) 



R — CO (Radiation condition) 

(3d) 


Note that the botindary condition 3a really expresses pressure continuity in 
the wake as well as no disturbance ahead of the wing in the wing plane so 
that 
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x> c— ytan 9 


(3e) 


di> 

at. 



= 0 


and 


<^=0 x<-ytane 


(3f) 


Inspection of both the flow tangency condition and the linearity properties' 
of the velocity potential permits further simplification of the problem. Let 
the equation of the surface be given by 

F(x,,X3,t) = X3-f(x,,t) = 0 ( 4 ) 


Then 


OF 

Dt 


at 


d^z 



(5) 


where 9«p/SX2 X2 -velocity induced by the wing and w is an imposed 
normal gust velocity. If the wing performs no unsteady motions, then the 
disturbance velocity potential must satisfy 


d4> 

dX2 


= u|^ +Wg{x,,X3.t) 


(6) 


Since the velocity potential is a linear function, superposition applies and 
a steady and unsteady component can be defined. 


<^ = <^3+^„S 


(7) 


Thus the boundary condition can also be separated into steady ^d unsteady 
relationships . 


dx2 <3X| 


d 4‘us 
dX2 


= WqCXpX^.t) 


(8a) 

(8b) 


This allows treatment of camber and angle of attack as separate solutions to 
the steady flow problem, while the response of the wing to an unsteady gust is 
obtained as simply the response of a flat plate to the gust. 
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Since the unsteadiness is represented by a gust harmonic in time, both 
the disturbance Telocity .potential and airload should be harmonic in time 


<^>{x,t)=<^{x) 

(9a) 

P{x,t) = P(x)e“^^ 

( 9 b) 

ACp(x,t) = ACp(x)e“*'* 

(9c) 


This boundary value problem can be converted to a singular integral eq.ua- 
tion which relates the local flow angle at the surface to the loading of the 
surface in a manner that expresses the flow tangency condition explicitly with 
the remaining conditions implied 


where 


w ACpie.T?) d^dTJ 

u ^ Jq 


AC p 


P(X2 = 0*)-P(X2=ffl 


( 10 ) 


and AC-jj "* ACp*, which is the result for an infinite wing in the sinusoidal 

7|-» eo 

gust field. A pressure coefficient increment due to the presence of a finite 
tip can he defined by 

ACp^ = ACp^-^P ( 11 ) 


which satisfies 


lim 

V-^oo ACpj.^0 


( 12 ) 


and -Which permits. Eq. (10) to he rewritten in the form 

cO C 

877-^e-'<M*V>=^ ^ K(x-e,y-^)[^Cp^(^,^)_^Cp^ (C,7?)]d^d77 


(13) 


Eq, (13) for points (x,y) on the wing can then be rewritten in the following 
form 
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( 14 ) 


87T^e-'(''xX-kyy) = ACp^.dCd-rj-J° f^KACp 
-CO 0 -oo'^o 

rODrC ^ ^ 

- /- / K ACp d^ dij 

•b *'0 

Now the upwasii induced by an infinite wing in the sinusoidal gust field is 
given by 

f 

f^f^KKCp dT 7 = 87 T ■ 

ico-'o ■ ^co ' u 

and the upwash induced on the setai-infinite wing by another semi-infinite wing 
extending from the tip to .T| = -« is given by the term 

/•orC _ 

iJo 

The two wings summed would make an infinite wing. Hence the integral eq^uation 
can be expressed as 

(e.‘^)dCci7?=-^^K(x-C,y-?)^Pco(e,^)dCd77 

y>o 

This equation expresses the physical phenomenon • that the presence of a finite 
tip induces an upwash increment over the wing surface' to account for the 
removal of that one half .of ah infinite wing which extends to The pre- 

vious result is a general result, applying equally to steady or unsteady 
motions as well as flat plates and cambered plates, and was also noted by 
Chu and Widnall (Ref. 8). 


Solution For The Pressure Field At The Tip Of The Semi-Infinite Wing 

The solution to the problem of determining the -loading on the semi- 
infinite wing makes use of the simplification afforded by the linearization of 
the problem. The angle of attack and camber contributions to the loading can 
by computed separately utilizing steady state< aerodynamic and then superimposed 
onto the unsteady results .of a- flat, plate responding to a sinusoidal gust. 

Equation 15 represents the basic equation to be solved for the finite 
tip pressure coefficient increment ACp^. The kernel function, K, has been 
calculated in a form suitable for numerical computation in Ref. 7 while 
ACp^ represents the solution to the infinite wing responding to the sinusoidal 
gust field and can be obtained from other sources (e.g. , Refs. 5 and 6). 
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One approach to a solution is a comhination of analytic and numerical 
techniques in which the known ACp^ is used to compute the upwash on a 
selected set of control points on the wing numerically. This is the right 
hand side of Eq. 15. Then the unknown ACp^ is computed at load points on the 
wing hy assuming a series representation made up of chordwise and spanwise 
loading fractions which satisfy the edge behavior of the semi-infinite wing 
.using 

yjTg X ytanfi (leading edge) 

^p X — C- ytanS (trailing edge) 

ACp — ^/v y-^ 0 

— ACp^ y — 00 

(where 6 is the distance to the edge of the wing), and a set of irnknown 
weighting functions. The weighting functions .are determined by satisfying 
Eq. 15 at the selected load points for which the right hand side was calcu- 
lated. This is essentially the method proposed in Ref. 8 and is quite satis- 
factory for small values of KM. As kM values increase, however, the chord- 
wise distribution of AC ' ' which defines the chordwise loading function, does 

pco 

not maintain as simple a form. 

The present solution begins by casting Eq. 15 iiito a numerical form 


0 I CO N — 

1 Z ACp^.. K(x-Cj.y-'^?j)AeiA^j = 1 ^ £ ACp^,^^K(X-e„.y- 77 ^)AC^A, 7 ^ ( 16 ) 

j=-oo i=i * rn=o n=i 


Thus, an infinite set of simultaneous equations, having a infinite number o'f 
unknowns can, in principle, be solved by standard matrix techniques for ACp , 
However, this equation can be ret/ritten as 

0 I . _ M N _ __ 

Z Z ^kli.^'^PcOi = ^^PcOnJ^klnm 

j=-co 1=1 ‘ (17) 

m=M n = l 

where M represents a practical limit on the spanwise position at which two- 
dimensional flow is attained, although ideally this is not satisfied except 
in the limit 

lim ACp- ,^ — 0 
M— CO 
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Thus- for Xa 2: L -* m & Mj AC-r, — AC and Bq.. 15 can be written 
’ 3 ’ Pnm P®nm 


01. _ MN OON 

’I Z Akii ACp +X I *klnm^Cp + ^ Z A^inm 

j=_cx)i=l m=0n=l M n=i 

CO N ^ 

i^^Z AkinmACpoij 

M n=l 


’’cOnm 


(18) 


M N 

= Z Z A 
m=0 n=l 


klnm 


^nm 


where now the whole left hand side represents the upwash supplied externally 
and the second terra on the right represents the effect of the two-dimensional 
inboard loading on the tip region. The entire region of the tip effect is 
hounded spanwise by m = o and m = M. Then Eg,. (l 8 ) may be rewritten in the 
form 


oON 

L.H.S.-IX 


Mn=i 


„ M N 

AKlnm^*^PcOnm”^ Z Akinm 

m=0 n=l 


ACp 


nm 


(19) 


which gives a set of N by M simultaneous equations which can be solved directly 
for the tip loading by standard methods. 

How, the remaining term on the left hand side of Eq. (I 9 ) may be written 
as 


CO N 

^ ^klnm 
M n=l 


P N 


AC 


^COnm 


= zz 

M n=I 


CO N _ 

Aklnm^CPcOnm"^^ ^_l^klnm^^PcOnm 


( 20 ) 


The term -\pnrn kernel function of subsonic aerodjrnamies (Ref. 7). 

BTormally the complexity of this function precludes analytic treatment, but 
for load points that are far removed from the downwash points, the kernel func- 
tion can he simplified so that some integration can he performed. Hence, in 
anticipation of these simplifications, the last term of Eq. (20) may be 
rewritten in integral form 


PN 1-00^0 

ZZAk,nm'ACp =Z Z ^klnm ACp^^^+gL J J K(x,-C, y,-^} ACp^d£d^ (21) 

M n=l """ M n=l -yp t) 

and Bq. (20), rewritten in the form of Eq. (21), can be determined partly by 
numerical and partly by analytic means. 

The formulation of the problem in this manner, with the external specifi- 
cation of the length, 1 , is a key development and permits a practical 
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generalization to the prohlem. Altering this length can alter the spanwise 
gradient In the tip region, although as this length increases, the effect on 
the tip loading will decrease. The choice of L for the present analysis is 
hased on the results of Eef . 8, and was taken to he 4 chordlengths for the 
calculations performed herein. 


The Pressure loadings On An Infinite Wing Yawed To A Compressible 
Stream Containing An hnbedded Sinusoidal Gust 

The solution to an infinite flat plate wing yawed to a compressible free 
stream containing an imbedded sinusoidal gast has been developed by Adamczyk 
{Refs. 5> 6). This solution represents an exact solution to the problem and 
was developed by casting the acoustic eq.uation as a modified Mathieu eq.uati’on 
and expressing the result, as determined by the boundary conditions, in an 
infinite series of Mathieu functions (Eef. 5)* The analytic form of these 
functions is .sufficiently complex that further analysis with then is difficult, 
but they are well suited to numerical computation, Ror ease in interpreting 
the equations of Ref. 6 the gust encounter geometry of Eef. 6 is given in 
Fig. k. 

The response of the infinite wing to a sinusoidal gust is characterized 
in Eef. 6 as a two-parameter solution. 

R- i^bCOSg 

|-M^COS^0 (22a) 


V sing losing r,MCos^v2_' 

I~m2cos 20 ' ■ 1-M2C0S26 Sina ' . 


kb 

] — M^COS^g 


cos^^ - sin2aj 


1/2 


(22b) 


The quantity, 3, is simply the chordwise reduced frequency component and has 
its direct counterpart in xmswept infinite wing theory. The parameter, Y, 
determines the behavior of solutions at large distances from the airfoil and 
is defined as the cutoff parameter in Ref. 6, It enters the problem formula- 
tion through a scaling procedure which reduces the acoustic equation to the 
following form (see Eef. 6) . 


*X,X,+ ^X2X2+ 


■$=0 


( 23 ) 


The fundamental solution to this equation is a source with an asymptotic limit 
given by 
14 



girR- 


R= VX|^+X2^ 


(24) 




Form its definition in Eq.. (22b), is either positive definite or negative 
definite, and hence Y will either be a pure real or a pure imaginary quantity. 
Now, i-f Y is imaginary the solution decays exponentially while if,Y- is real 
the solution approaches the asymptotic form of a cylindrical acoustical wave 
propagating outward from the origin. A key parameter here is M*, the phase' 
Mach number of the disturbance along the span relative to the freestream 
flow (see Fig. 5a). Whether or not the cut-off parameter, y, is real or 
imaginary depends on whether, or not M*- is greater or less than one. 

Kie effect of the relative phase velocity, M*-, is to introduce a spanwise 
velocity for which the flow appears stationary in an axis system also moving 
at the spanwise velocity. The nodal lines of the. gust in this moving axis ' . 

system are also stationary. The problem is thereby transformed from the non- 
stationary problem of a gust encounter at some angle relative to the vnlng,, to 
a problem involving a stationary gust front parallel to M* and lying across 
a wing idiieh moves spanwise at a velocity M* cosa. 

Adamczyk (Ref. 6) has successfully simplified the exact solutions for 
distxirbance pressure for limiting cases of this compressibility cut-off param- 
eter. From Ref. 6, the disturbance pressure can be written in several differ- 
ent forms for various values of the cut-off paramter, Y. When Y^ is a small 
positive quantity, 


Ps = - 


/>qUu 2 cosg 

(l- ^ i + X| 


l+Z 


•] 


(25) 


exp{i[-/3M2cos2ex,+ X3-kUcos(0-a)t]} 


with 


iHo(3) + i (y/ffl^[Ho^'ky/2)-Hot» (ff). 


,(l)i 


0< y2<0.l 


When Y^ is negative, 

p . po UO2COS e / l=^ H(y,jg) 

(1-m2cos20)'/2 V |+X| ro(y)+l,(y) ‘ 


exp[-(y +\pM^cos^e^ + (|_m5^cos2 0 )i /2 X 3 -ikUcos( 0 -a)t 


(26) 
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where eyipl-\-y/y^+ y9^ cosS S)(l+^^^j/( !+ 2 7 rV^y 24 -^ 2 (H- ^ )] 

H(r.^)= ^ 


\+ir^yZ^j^Z [l + C0S^S+ir-/y2 + ^2 sinS] 

8 = tan-'(y//3) 0>Y^ 


For larger positive and negative values of Y^j 


P /?>UU2COS0 

's 




where 


(,_M2cos 25)I/2 I yTxr 2 

exp-[i[(y-y 9 M^cos^ 0)x(^+l^b sin a X3-ui<cosC0-a)tj| 

x-;=x,+i |y^|>0.7 

E(x)=C(x)-iS(x) 


(27) 


e 1 x)= CONJUGATE E(x) 

The corresponding pressvire loading on the yawed infinite flat plate wing 
responding to the sinusoidal gust can be obtained from 


r 2Pq - 

ACp^=CONJ|^ — j — expi Ukcos(0-a)t 


^CD 


JL p 

2'° 


(28) 


The complex conjugate is taken since Adamczyk solved the conjugate problem to 
the formulation described herein. 

In the present analysis Eq. (25), (26), or (27) will be used, depending 
on the value of y^. If Y^ does not fall within the prescribed ranges, the 
exact solution will be used (cf . , Eef. 6). , * 

The steady flow contribution from a flat plate at an angle of attack is 
obtained as a limiting case of the unsteady analysis at zero reduced frequency. 
The effect of the steady flow camber problem can be treated by obtaining a 
potential flow solution for an infinite wing and then accounting for compress- 
ibility by way of the Karman-Tsien relation (Bef. 9). A suitable potential 
flow result for an arbitrary mean line is given by the Munk integral (Ref. 3). 
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where = cos“^ . Instead of using the Prandtl-Glauert correction factor 
to account for compressibility effects, the Kaxman-Tsien relationship, which 
is more accurate at high' subsonic Mach numbers, is used. This relation is 
given hy 


CpM ! 

Cp, yi-M2 + Cpi 

l+yi-M2- 2 


(30) 


where C is the pressure coefficient at the Mach nimber M, and C is the 
incomipressihle pressvire coefficient at the same location in the flow. 


At this point attention can be turned momentarily to the problem of 
determining the pressure in the flow field near the wing but off the sxirface. 
■The leading edge pressure singularity, which is inherent in the theories of 
infinitely thin wings and which generally causes no problems in computing the 
aerodynamics of these wings, must now be treated. If the loading on the wing 
is represented by pressure dipoles, the disturbance pressures at any point in 
the flow may be directly determined by 


1 a 


r - M{-x-e) + y{x-C)^ + (l-M2)[:x2+(y-77)2j -j 

4v dXp R 

WING AHEA 



where / I I I 

(x-0 + ^ 2 ^ +{y-^) 


(31) 


Unfortunately, the infinite strength dipole at the leading edge will then 
present infinite pressures throughout the fluid, an obvious physical unreality. 
Hence, the loadings in the vicinity of the leading edge must he modified to 
eliminate this singularity. 


Lighthill (Ref, 10) has successfully modified incompressible two-dimen- 
sional thin airflow theory to eliminate the leading edge singularity. It is 
required that the flow near a smooth leading edge be the flow over a parabola, 
while at points away from the leading edge the flow is adequately predicted 
by thin airfoil theory. This moves the singularity inside the airfoil contour 
and provides a smooth flow exterior to the thin wing. As a result the pres- 
sure loading on a thin wing can be given, to a first approximation, by the 
following equation, 

ACp _ X 

° (32) 


where x is the distance from the leading edge, AC^ is the actual pressure 
loading across the wing, AC ^ ^ is the thin airfoil theory result at the 
corresponding chordwide position and p- is the airfoil leading edge radius of 
curvature. Since ACp,p ^ varies as (x)“V^ near the leading edge it is easily 
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seen that ACp approaches zero continuously as the leading edge is neared. As 
X approaches unityj approaches its thin airfoil ‘value. Without direct 
proof, it is assumed the same results will hold for subsonic compressible 
flow; this is based on the fact that the behavior of the leading edge singular- 
ity is the same regardless of cotapressihility or unsteadiness. That is, for 
■unsteady compressible flow ACp varies as near the leading edge. 

The solution to the problem of predicting the pressure field in the 
vicinity of a swept, semi-infinite wing consists of a numerical solution to- 
Eg.. ( 19 ), utilizing the proper value of AC . The necessary dot* let-lattice 
influence coefficients were developed by Casing, et al (Ref. 7 ) and the 
underlying assumptions in their derivation restrict the solution to the case of 
tip chord parallel to the flow. Once the thin wing results for the loading 
are determined, the leading edge singularity is removed by Eq. (32) and the 
pressures' in the near .field are determined from the dipole distribution, Eq-. 

( 31 ). ' - ' ■ ’ - ■ - 


An Asymptotic Solution For The I^pwash Induced At The Tip By 
Portions of' The Wing Far Removed From The Tip 

Equation (21) which is repeated below for convenience, describes the 
upwash on the tip due to that portion of the wing which experiences two-dimen. 
sional flow. 

E 1 Akimn ^ E Akimn ^^PcOmn + 

M n=l M n=l X. -b 

Here A^-| represents a numerical form of the- unsteady subsonic flow kernel 
function, x„- 5.), both of which are given in Ref. 7- The continuous 

form, Kj of the kernel function is more amenable to analytic treatment. The 
following relations, from' -Ref. Tj are required to perform the integration in- 
Eq. (21). 

K(Xo,r) K|/r^ K, = I,+ — 

R.Cx„2+/S^yo^)''^ 

Bie 'inte^ai; is given exactly by 


r- Yo- 


Yo =y- 77 


(33) 
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, r” e-"’i"du 

(l+u2)V2 (34) 


and approximately "by ' 


where 




U| 

(l+U,2)'/2 


ik|Io 


(Ui.ki)] 


(35) 


Io(U|,k|) 



^^g-ncu| (nc-ik|) 
ri2c^+k|^ 


(36) 


was developed by Laschka (Eef. 11) and is reproduced in Ref. 7. !Phe Lasehka 
coefficients necessary for this approximation .are n, and c. 


In general the spanwise variation in AC is oscillatory due to the span- 
wise cotoponent of the gust. From the -form or Eq.s, (25-), (26)-, or (27) it c.an’ 
be seen that 


^Cpj{£,) e' 


I'H' 


iySM2cos20^i*ik sinnCj 


(37) 


Thus the term in integral form' on the right hand side of Bq.. (21) can be 
written as 


r ^ ^ g-i/SM^ cos 20 ^, r^^giKsinaiCa g- ^i>‘rCl>cos 0 + (xj-^^jsine] 

4 ■' ■ 4 



Defining 




this can be restated as'. . -- 


— /^icp g ifcsinaxj ^-^(x,-^,)cos0 'tks!nd*QSin0Ji7r.t^,j^ 

Define the distance, along the span between the point at .which the upwash is 
desired and the load is applied- to be 11^., If this .quantity is large certain 
approximations in the kernel function, K^/r^ can be made. In, particular, 


■'r^ = yo^ = ‘^r^cos^0 
rR = i?r2cos2e{tafi2(-.0-i-€) 

M(tang (--0 +g) +/3^)'^^_ ton (-04-^) ■ 
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Combining these approximate forms, the kernal function has the following 
approximate form 


iSl 1 

r2 Lc4 CiJ 


g-iC2% , N bne'*^2'^r 


+ J_ y .. r:* . : — i. y — 

n. ^ -J O C4 ’Jr 


’Jr' 


C- n.| 


I ^ _dne 

_ m 

4 n=i 


-iC2’7r 


with = cos^6 
Ui 


Co - TT cos 9 


C3=l- 


Vi + U|^ 


'2 U 
C4'C| 


[gn^+’Jr^] 
[M(tan^ {- 0+^) + /3^)''^^— tan(-0 + 6)j 


(38) 




bn=Qne 


dfi“ bfj * 9n 

U| s (tan^ (- 0 +e ) + — tan (~ 0 + e )j /^2 


9n= 


nCU 

U>COS0 


The angle, e, is the angle between the load point on the wing and the 
point at Which- flow tangency is satisfied, measured from a spanwise line 
through the load point (see Fig, 5h). Note that this angle is a function of 
chord, but for large spanwise distances between the loadpoint and the flow 
tangency point this "variation can be neglected. The integral 

f ® iii dv 

•'x,-X 

can be evaluated by first defining a parameter 

Oisiii0 WCOS0 r M(tan^(-0+e) + /3^)''^^— tan(-0+e) i 


D=ksina + 


U 


U 


P‘ 


Use of this parameter then permits evaluation of the following integral forms 
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-00 e 


-iD77r 


Xj-X 
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,-co p'-iO’Jr 


-00 
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e" 'D’Jr 


>^3-^ + 




The integrations are performed in the complex plane and the following simple 
form is obtained 


/, 


rOO 


-iD’Jr K| 


X3-X 


<1 d-r, 1 

Fa ‘‘’>' = Lc4'^c7J 


One 


l2 SGN(D) it ^ 

2 2 n=, Qn 


,-ncu, 


j^l- SGN(D).2e 


-iGnI 


1 


(39) 


where 


Gn = Dg-n 


This calculation represents a small, real component to the upwash due to 
portions of the wing far removed from the tip region. 
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MJMERICAL RESULTS FOR INVISCID PRESSURE FEEIL 


The computational procedure to determine the tip loading consists of 
representing the tip region (0 L) "by the doublet -lattice model des- 

cribed in Ref. 7 . Use of this model presently restricts the tip chord line 
to be parallel to the freestream. The length, L, at which the tip influence 
substantially disappears is adeq,uately given by chord lengths. The inboard 
reagion of the wing is described by a combination of analytical and numerical 
forms, utilizing the infinite wing results of Ref. 6, the numerical kernel 
function of Ref, 7 , and Eq., ( 39 ). Camber effects due to an infinite wing in 
steady subsonic flow are treated by combining the incompressible solution 
of Eq. (29) with the Karman-Tsien equation (Eq. 30) to account for compres- 
sibility. Equation (I9) is then solved numerically for the unknown pressure 

loads on the tip, AC . The final loads on the tip region are obtained by 
^ ^ pnd / _ \ 

removing the leading edge singularity from the ACp^j^ via Eq. ( 32 ). 


Before applying the approximate pressure fields described by Eqs. ( 25 ), 
(26)5 and (27) to the solution for the semi-infinite wing, it is desirable to 
establish some limits of accuracy for these analytic forms. The -quantity of 
primary interest is, of co\irse, the pressure. However, the lift transfer 
function Tj^ (which is a normalized lift coefficient that includes the effects 
of compressibility and the wake) is also of interest and has been examined 
in the figures which follow. Adamezyk (Ref. 6) has developed these transfer 
functions for the following tentative values of the cut-off parameter. 

T^_^ibkC0Sa^ Jo(/3)-iZJ|(g ) ][ COs20)+iJ|(M2;SCOS^0}] 

- i>r>'o ( 40 ) 


T| = 


j^(r+i^M^cos^0)+Ii(X + 9) ^^giUbcosa 


io(r)+^i(/) 


[|-s][l + ^sinS 


H(r,/ 3 ) 


|cos 8 - l+2^vyr2+/32 Cl+a 5 Sins ' ll 

l + vy/ 2+/32 [l + cos^ S +vy 72+^2 sinS] 


y2<o 


(Ul) 


1 




Q 2 i[y-l 3 M^cos^ 9 ) 


[y2(y-;QM2cos2 0)] ■ 


t[y^] 


2(y— COS^ 0) y-^IVl^COS^0 


g2l(y- ^M^COs20 ) 

y-^M^cos^0 yy + y3M2c^e 


i y 2 y 
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where Z and H (YjP) are defined in Bq.s, (25) and (26). Note these transfer 
functions are referred to the gust amplitude at the leading edge rather than 
mid- chord . 

In Fig. 6 the "exact" solution for is compared with the approximate 
forms for positive in Eqs. (40) and (42). (The exact calculations were 
obtained hy a direct numerical integration of Eq. 21 of Ref. 6.) The upper 
panel contains the comparison for M = 0.3 siid the lower panel for M = 0.6, 
both over the range of reduced frequencies 0.05 ^ kb s 6.0. In both cases it 
is seen that Eq, (4o) which is the low frequency approximation, compares 
favorably with the exact solution up to kb = 0.2. Beyond this value Eq. (40) 
diverges from the exact solution, with a greater divergence at M = 0.6 than at 
M = 0.3. In contrast to this, Eq_^ (42), which is the high frequency approxi- 
mation, is in poor agreement for_kb <0.5, but improves in its ability to • 
represent the exact solution as kb increases beyond this value. This is also 
illustrated in Fig. 7 iii which the magnitude and phase angle of T^^ are plotted 
versus the cut-off parameter Y. Based on the results of these two figures, it 
appears that Eq, (4o) can he used for Y < while Eq, (42) can he used for 
Y > 0.3. 

Figure 6 shows the chordwise pressure distribution on an infinite wing 
at a Mach number of 0.3 for a reduced frequency of 0.1 in part a (and shows 
the variation with reduced frequency at x/c = 0'.2 in part b). Comparison 
of these results indicate that Eq. (25) yields acceptable pressure distribu- ' 
tions for y < 0.2 while Eq, (27) is more than adequate for Y > 0,2, Figure 9 
shows a similar set of comparisons for a M6ich number of 0.6_^ For this set of 
conditions, Eq. (25) re£resents an adequate solution up to kb = 0.1 and is 
possibly extendable to kb as 0.2. Equation (27), however, appears to compare 
quite well with the exact results beyond kb — o.2. This implies that Eq. (27) 
still gives an adequate descritpion of the pressures for Y > 0.2 while Eq. (25) 
is valid for Y < 0.2. 

It should be noted that the foregoing comparisons relate to conditions 
for Y^ > 0, in which acoustic propagation is present. The conditions for 
Y^ < 0 have not been investigated herein and remain a field for further study. 

Deviations begin to appear in both the lift transfer function and the 
pressure distributions at high values of reduced frequency at both Mach numbers 
considered with the deviations becoming more pronounced as Mach number 
increases. This is a result believed due to the truncation of the infinite 
series of Mathieu function which describes the exact solution. Equation (27) 
is a solution comprised of the first two terms of an infinite series in reduced 
frequency in which the disturbance due to the plate is modeled by acoustical 
waves emanating from the leading and trailing edges. This solution inherently 
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increases in accuracy as the reduced frequency increases (at a given Mach 
number ^ 0) and so it is expected the results of Eq. (27) are more accurate 
than the truncated representation of the exact solution at high values of kb.. 


Pressure Loadings on The Semi- Infinite Wing 

Figures 10 and 11 describe the spanwise variation in chordwise pressure 
distributions of .the semi-infinite ving at various sweep angles. The wing is 
an uncambered flat plate at an .angle of attack of O .25 rad in incompressible 
flow. Generally, the chordwise distributions are of similar form, and the 
spanwise variation decreases smoothly to zero at the tips. The effect- of sweep 
is to lower the overall magnitudes. 

Figure 12 represents the steady response of a cambered (WACA 4415 profile)- 
unswept semi-infinite wing at an angle of attack of 0.25 rad. A comparison, of 
this figure with Fig, 10 shows that the effect of camber is -to increase the 
£hordwise load distribution. For completeness, the unsteady distribution for 
kb = 0.1 is Included in Fig. 13. As stated earlier in the text, the complete 
solution Is the linear superposition of the steady and unsteady components of 
the response which are illustrated in Figs. 12 and I 3 . 


Pressure Distribution m The Vicinity of Tip of Semi-Infinite Wing 

Figure l4 shows the pressures in a region surrounding- the tip of a semi- 
infinite unswept wing in steady sxibsonic flow at a Mach number of 0.2. The 
angle of attack is taken to be 0.1 rad and camber is zero. In each panel the 
pressure coefficient is plotted horizontally versus the normal distance to the 
plane of the wing. The left column of panels at y/c = 2.0- is the distribution 
two chordlengths inboard of the wing tip, and the right column of panels at 
y/c = -2.0 is the distribution two chordlengths outboard of the tip. The 
first horizontal row at x/c = 0 .I 5 is the distribution at points that are I 5 
percent chord aft of the leading edge locus, the second- row for x/c = 0.95 
(consisting, of -a single panel)- is at a point that is 95 percent chord aft of 
the leading edge locus, and the third row at x/c = 1-45 is for points that 
are 45 percent chord aft of the trailing edge locus. The right hand column 
can he .dismissed immediately as representing typical variations of pressure 
in the vicinity of, but external to, a lifting surface and its wake. The top 
panel of the left h^d..column for a point near the wing leading edge contains 
the characteristic to the dipole loading. Negative pressures (Cp < 0) 

are observed f or' x^' > ’0 Vhjtle positive pressures are observed for Xg < 0. The 
bottom panel of the left hand column for a point in the wing wake shows the 
far field dipole influence at points far removed from the wake {x 2 |» 0 . 
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However, as the wake is approached a reversal of this trend is observed and 
the boundary condition of zero pressure difference in the wake is satisfied. 
Finally, the center panel of the left hand column represents a point near the 
wing trailing edge. Here the pressure variation is transitional between the 
dipole Jump on the wing and the continuous pressure in the wake. As the wing 
is approached the pressure magnitude initially increases. At a small but 
finite distance from the surface of the wing, the trend is reversed and the 
pressures tend toward zero, but terminate in a discontinuous dipole Jump on 
the surface. This combination of form, near a point where the boundary condi- 
tions change abruptly, is characteristic of the elliptic equations which 
describe subsonic steady flow. 
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VISCOUS FLOW CONSIDERATIONS 


General 

Although, an accurate solution of the potential flow eq.uations can give 
a qualitative picture of the flow in the region of an airfoil tip, the flow 
field in this region is affected significantly hy three-dimensional viscous 
phenomena arising from boundary layers on both the upper and. lower wing sur- 
face as well as on the wing tip itself. Viscous effects could he taodeled 
partially in a potential flow solution through use of isolated and discrete 
vortex filaments. The filaments would represent the vorticity originally 
generated in the viscous hoimdary layers which break away from the airfoil 
to form the wake. However, a complete potential flow solution which includes 
discrete vortices would require a calculation of the interaction of the 
vortices with each other as well as their combined effect upon the wing. 
.However, even such a complete potential solution could not satisfy the suxr- 
face no-slip boundary condition. Although failure to satisfy the no-slip 
boundary condition may not be serious in determining the pressiire field, it 
may lead to inaccurate predictions of secondary flows in the vicinity of the 
wing even in regions which have little direct viscous effects. In addition, 
a vortex filament solution would require assumptions on the vortex core size 
and upon the position and angle at which the filaments leave the wing as both 
of these items are dependent upon viscous effects. 

Since a complete interacting vortex filament-potential flow solution 
would not be able to model either the airfoil boundary layer or the viscous 
core correctly but would still present a formidable computatibnal problem, 
it is appropriate to consider viscous solutions of the tip flow field. One 
possible viscous solution would attack the full three-dimensional Wavler- 
Stokes equations. Solutions of the full Navier-Stokes equations can now be 
considered possible (e.g., see Briley and McDonald (Ref . 12)), however, their 
relatively long computer run times dictate that the full three-dimensional 
Navier-Stokes equations be solved only when no suitable alternative is avail- 
able. 


Although flow situations exist for which a three-dimensional Navier-Stokes 
solution is appropriate and necessary, the wing tip problem does contain sim- 
plifying features that may alleviate the need for a full three-dimensional 
Navier-Stokes solution. In particular, the tip flow region contaiiis a primary 
flow direction which, for the zero camber case considered in the present 
effort, is in the plane of the oncoming flow and at only a small angle rela- 
tive to the wing chord. Furthermore, the velocity component in the streamwise 
direction is considerably greater than velocity components normal to this 
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streamwise direction and derivatives in the strearawise direction are expected 
to be considerably less than derivatives normal to this strearawise direction. 

In this sense, the problem appears to be similar to the classical three- 
dimensional boundary layer problem. However, the classical three-dimensional 
boundary layer approach imposes two additional assumptions which are not 
valid in the airfoil tip region. Under the first of these additional assump- 
tions the static pressure is assumed to be invariant in a plane normal to the 
strearawise direction; obviously, such an assumption upon the static pressure 
clearly makes classical three-dimensional boundary layer theory invalid for the 
airfoil tip problem. Secondly, three-dimensional boundary layer theory assumes 
that significant gradients exist in one direction only, i.e., that a thin 
relatively flat boundary sheet exists wherein viscous effects resulting from 
the mean flow gradients are significant only as a result of the velocity 
change across the sheet and not along the sheet. Obviously the assumption 
of a thin relatively flat boundary sheet is inappropriate iri the region of 
a wing tip. 

The need to develop a three-dimensional forward marching viscous calcula- 
tion procedure which is more general than three-dimensional boundary layer 
theory has been motivated by a variety of fluid mechanics problems such as 
three-dimensional duct flow problems, three dimensional jet problems, and the 
airfoil tip problem. Therefore, considerable recent effort has been expended 
upon development of computational procedures for three-dimensional viscous 
<• flows with a dominant strearawise direction. These procedures would not be 
limited by the classical three-dimensional boundary layer assumption of con- 
stant static pressure in planes normal to the approximate flow direction but 
would still allow treatment of the flow field in question as an initial value 
problem. Such three-dimensional viscous forward marching solutions have been 
developed by Patankar and Spalding (Ref. 13), Caretto, Curr and Spalding 
(Ref. l4) and Brilej'’ (Ref. 15) all of whom developed numerical solutions for 
laminar incompressible flows in straight ducts with rectangular- cross sections. 
The governing equations were solved by integrating in a primary flow coordin- 
ate direction while retaining viscous stresses in both transverse coordinate 
directions as opposed to only one direction for three-dimensional boundary 
layer^heory. In addition, certain assung)tions were made about the behavior 
of p'ressure gradient terms for incompressible flow to permit solution by 
forward marehing integration. Subsequently, this general approach has been 
used to compute laminar incon 5 >ressible flow in helical tubes by Patankar, 

Pratap and Spalding (Ref. l6). 

Recently in companion studies, Briley and McDonald (Refs. 12, I 7 ) have 
developed a stable and efficient noniterative implicit numerical technique 
for application to systems of coupled nonlinear multidiinensional parabolic 
and/or hyperbolic equations. These general techniques were applied in 
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Ref- (17) to the computation by forward marching integration of laminar 
supersonic flow in rectangular jets. Finally, the technique was extended 
in Ref. (18) to the prediction of subsonic turbulent compressible flow in 
curved ducts. In the present effort this same basic viscous, three-dimen- 
sional forward inarching technique is applied to the airfoil tip problem, and 
a preliminary assessment of its potential for calculating the tip flow field 
is made. 


The Governing Equations 

Central to the present analysis is the formulation of approximate 
governing equations which can be solved by forward marching integration in 
the direction of a "primary flow". The entire flow field then can be 
obtained by a sequence of essentially two-dimensional calculations, and this 
feature of the method results in a substantial saving of computer time and 
storage compared to that which would be required for solution of the full 
(elliptic) Wavier-Stokes equations. Although the present effort utilizes a 
cartesian coordinate system to represent the airfoil tip flow field, the 
equations have been derived and coded in a more general orthogonal system. 

The reduction to a cartesian system is straightforward and is indicated 
during the discussion. 

The equations are derived in a curvilinear orthogonal coordinate system 
which is aligned with the flow geometry such that one coordinate direction can 
be identified as the primary flow direction while the remaining two coordinate 
directions determine the secondary flow plane. The transverse plane is 
assumed to be perpendicular to the airfoil. Only the tip region of the air- 
foil within one-half chord of the edge is considered. The airfoil tip is 
assumed to lie in the plane 2 = 0 between 0 ^ x £ c and -h/2 ^ y ^ h/2 where, 
h is the airfoil thickness . The leading and trailing. edges of the airfoil are 
assumed to lie in the planes x = 0 and x = c, respectively. A sketch of the 
airfoil tip coor^dinate system used in the viscous analysis is presented in 
Pig. 15. 

In the general case x, y, and z represent the approximate streamwise 
(primary flow) and two transverse coordinates in an orthogonal coordinate 
system, respectively. Since the analysis consideres a general orthogonal sys- 
tem, metric coefficients h^^, hp, and are defined such that an incremental 
distance (§s)^ = (h^^bx)^ + (h2by)^ + (^3^2)^; for a cartesian system hj_ = hp = 
= 1- The governing equations are derived from the Navier-Stokes equations 
describing the compressible flow of a viscous, heat conducting, perfect gas. 

In vector form, these equations are given by 
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a^/at + v . ipu) =0 


(J^3a) 

y 3 au/at+/?(u- v)u +vp =F (^3b) 


where U is the velocity vector, p is density, t is time, p is piessure, F is a 
vector representing the viscous force acting on the surface of a fluid element 
and V is the' gradient operator . The total enthalpy E may be specified through 
an energy conservation equation (Ref. l8); however, in the present effort the 
assumption is made that the total temperature is constant throughout the flow 
field. The eqioation of state is p = pRT, where R is the gas constant. Eiipres- 
sions defining F in an orthogonal coordinate system are given by ’Pai ('Ref. 19) 
and may be found in Ref, l8. 

In the present analysis the reasoning followed to estimate the order of 
magnitude of viscous terms and turbulence quantities is the same as is often 
employed to derive three-dimensional boundary layer equations.- The rationale 
is discussed in detail in Ref. l8, and fcr convenience this discussion is now 
repeated in a condensed form. As discussed in Ref. l8, it is assumed that, 
viscous effects are negligible except in thin layers near solid boundaries, 
and thus boundary layer concepts can be employed to examine the relative 
importance of viscous terms in the governing equations. Consequently, vis- 
cous terms which are considered important for boundary layer flow on solid 
surfaces aligned with either of the two coordinate planes which are parallel 
to the primary flow coordinate are retained throughout the entire region of 
the viscous flow calculation; -other viscous terms are neglected. In the 
context of the present investigation the lower and upper surfaces of the wing 
represent solid surfaces aligned in the ’y '-'^.irection and the tip represents 
a solid surface aligned in the ' z '-direction. Both of these solid surfaces 
are parallel to the primary flow direction, 'x'. If the flow is turbulent, 
the governing equations are time -averaged in the usual manner for' turbulent 
flows (e.g., Hinze, Ref. 20). The dependent variables are represented a's the 
sum of a time-averaged quantity denoted by an overbar ( ) and an instanta- 
neous fluctuating quantity denoted by a prime (’}. This process of averaging 
produces turbulent correlations which are conventionally termed Reynolds 
stresses. The order of magnitude of viscous terms, including the turbulent 
Reynolds stresses, is examined under two sets of- circumstances, name'ly, those 
appropriate for boundary layer flow near either of two types of wall, one for 
which y is constant (a y- wall) and one for which z is constant (a z-.wall). 
Hear a y-wall, (v/u^Ep) and 5 ( )/5(y/%EF) are assumed to be of order '6 and l/6 
respectively; near a z-wali, (w/U^gp) and 5( )/3(2 /Lj^f) are of order <5 and 
l/5 . Here , 6 is the shear layer thickness. If the viscous terms are to be 
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of tlie same order of magnitude as the remaining terms, then the dimensionless 
molecular viscosity must be of the order 6^. Similarly, it 

is assumed that dimensionless turbulent double correlations are of order 6 
(i«e., turbulent fluctuations are of order 62 ) and that triple correlations 
are negligible. All other dependent variables and derivatives are assumed 
to be of order unity. Retaining only those terms which are of order unity 
in either of the two types of shear layers described above, the viscous terms 
are greatly simplified, and the time -averaged governing equations can be 
written in the following form: 
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For entirely supersonic flows, Eq.s. (4Ua~d) together with boundary and 
initial conditions can be solved without further apprcKimation by forward 
marching integration in the x direction as was demonstrated by McDonald and 
Briley (Ref, 17) for laminar flow in rectangular jets. For siR>sonie flow, 
however, the inviscid flow region is known to be governed by equations which 
are elliptic; that is, by equations which require downstream, boundary condi- 
tions for solution. Therefore, in the subsonic case, if the pressure is 
required to emerge from the solution it is by no means clear that a stable 
numerical solution could be obtained from a forward marching calculation and 
even if a stable numerical solution could be obtained, solution by forward 
marching integration is not physically appropriate, at least not without some 
sort of iterative procedure to satisfy the downstream boundary conditions. 
Therefore, it appears that in the subsonic case some approximation must be 
made in regard to the pressure field. The present airfoil tip analysis follows 
the approach of Ref, l8 and assumes the pressure field appropriate for invis- 
cid flow represents a given reasonable first approximation to the actual pres- 
sure field. . Thus, the inviscid streamwise pressure gradient computed with 
appropriate downstream boundary conditions is "inqiosed" upon the flow, as a 
known source term in the strearawise mojnenttim equation much as in conventional 
boundary layer theory, so as to permit solution by forward marching integra- 
tion of the viscous flow equations for subsonic flows. This choice of the 
pressure approximation is not the only possible choice and the manner in which 
the necessary pressure approximation is made is still a subject of current 
investigation- In any case the imposition of inviscid pressure gradients 
incorporates a priori the elliptic effects associated with a subsonic pres- 
sure field without the necessity of solving elliptic equations other than for 
an inviscid flow. In the present application, the inviscid pressiare distri- 
bution was obtained via the inviscid analysis described in the previous 
sections . 


Method of Solution 

The governing equations can be solved (after modeling the Reynolds 
stresses in the case of turbulent flow) following the general approach 
developed by McDonald and Briley (Ref. I 7 ) for laminar supersonic flow in 
rectangular jets. A detailed discussion of the calculation procedure is not 
included here, as such a discussion would be lengthy, and discussions of the 
general approach are available elsewhere (Refs. 12, I 7 , I 8 ). The method used 
is based on an implicit scheme which is potentially stable for large step 
sizes. Thus, as a practical matter, stability restrictions which limit the 
strearawise step size relative to the transverse mesh spacing and which become 
prohibitive for even locally refined meshes (e.g., in laminar sublayers) are 
not a factor in making the calculations. The general approach is to enqploy 
an iinplicit difference formulation and to linearize the implicit equations by 
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expansion aPout the solution at the most recent streamwise location. Terms 
in the difference equations are then grouped by coordinate direction and one 
of the available alternating-direction implicit (ADl) or splitting techniques 
is used to reduce the multidimensional difference equations to a sequence of 
one-dimensional equations. These linear one -dimensional difference equations 
can be written in block-tridiagonal or a closely related matrix form and 
solved efficiently and without iteration by standard block elimination 
techniques. The general solution procedure is quite flexible in matters of 
detail such as the type and order of accuracy of the difference approximations 
and the particular scheme for splitting multidimensional difference approxi- 
mations. Based on previous experience, however, it is believed that the con- 
sistent use of a formal linearization procedure, which incidentally requires 
the solution of coupled difference equations in most instances, is a major 
factor in realizing the potential favorable stability properties generally 
attributed to implicit difference schemes. 

As indicated above the details of the method of solution are quite 
lengthy and have been discussed In elsewhere. Rather than repeat the discus- 
sion in detail within the body of the present report, the specific numerical 
method is described in Ref* (l8) and given briefly in Appendix A. The topics 
in this appendix include the difference operators, the linearization process, 
the difference equations and the matrix inversion procedure. 

In brief, the numerical technique represents first and second derivatives 
by either first or second order difference formulas. However, a single 
direct substitution of the difference representations into the differential 
equations would lead to a nonlinear set of algebraic difference equations. 
Therefore, a linearization procedure is required. The linearization procedure 
is based upon previous work by ‘Briley and McDonald (Refs. 12 and 17) which 
assume the solution to the equations is Taylor expandable and, -therefore, 
approximates the solution at streamwise station n'+ 1 through a Taylor expan- 
sion of the known solution at station n. The result is a set of linear 
difference equations representing the nonlinear differential equations at 
station n + 1. The finite difference solution itself is an ADI solution 
based upon the splitting technique of Yanenko (Ref, 21) which results in a 
two step procedure in advancing the solution from station n to station n + 1. 
Each step requires the inversion of a block diagonal matrix which is accom- 
plished through a standard block elimination technique (e.g,, Ref. 22), 
Eurther details of these procedures are found in. Appendix A. 
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Grid Transformation, Difference Representation 
and Boundary Conditions 


The minimization of the computational effort necessary to compute a 
solution requires that a nonuniforra grid spacing he used to ensure that grid 
points are closely spaced in regions "where the solution varies rapidly and 
widely spaced elsewhere. In the present calculations, steep gradients occur 
in the boundary layers on the airfoil surface and in the vicinity of the 
airfoil tip and it is in these regions that closely spaced grid points are 
required. In the present procedure grid point packing is obtained using a 
transformation originally devised by Roberts (Ref. 23). For the purpose of 
demonstrating the Robert's transformation suppose that N grid points are to 
bemused in the range y^ ^ y ^ and a boundary layer or sublayer of thick- 
ness ®(y„~y. ) is present near y, , then Roberts' transformation tl(y) is given 
by 

T(y).N + (N-l) log (.ag.) ^^5, 

where a = ^ = a /(l"e)j and c = y^. The use of equally spaced points 

in the transformed coordinate T] ensures. an adequate resolution of both the 
overall region between and y and the boundary layer region between y^^ and 
y, + ® (yo“y^)‘ Derivatives with respect to the physical coordinate y are 
obtained from the following formulae; 


'^-[171 
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The use of three-point difference operators for 'll derivatives in Eqs . (46) 
produces similar operators for y derivatives. These y-derivative operators 
can be computed at the start of the calculation and stored, along with the y 
locations of grid points. 

In the present calculations the grid was resolved in the y-direction 
along both the upper and lower surfaces of the airfoil and in the z-direction 
in the vicinity of the airfoil tip. In the y-direction the limits of the 
computational region were taken between -.5 ^ y/c ^ .5 and transformations 
were performed both above and below y - 0 (the wing eenterplane) . For the 
transformation above the wing eenterplane, y^ was taken as 0.5c, yj^ was taken 
as 0.0 and e was taken as .04; below the wing centerline y^ "was taken as -0.5c 
yj^ was taken as 0.0 and e was taken as .04. Similarly, in the z-direction 
transformations were performed both inboard and outboard of the wing tip plane 
In these transformations y^ was set to zero, e to 0.04 and y^ was set to -fO.4 
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Re$ults 


Inviscid Calculation - Although the present section is concerned 
primarily 'With results of the viscous calculation procedure, the viscous 
procedure requires as input an inviscid approximation to the pressure field. 
Since this inviscid pressure field is used as an imposed source term in the 
streamwise momentum equation, the inviscid pressure distribution plays a 
significant role in determining the viscous results. Therefore, before 
examining the results of the viscous procedure, it is appropriate to review 
briefly the inviscid calculation procedure and to examine the inviscid pre- 
dictions in the airfoil tip region. 

Under the present effort, the inviscid flow field was determined by a 
doublet lattice method (see previous sections of the present report). For 
steady flows the doublet lattice method reduces to the modeling of the wing 
by a group of horseshoe vortices. The wing is divided into FI -rectangular 
sections and the bound filament qf a horseshoe vortex is placed in each seg- 
ment. In addition, one collocation point is placed in each rectangular seg- 
ment. Generally, the bound filament is placed along the local quarter chord 
line of each rectangular element and the collocation point at the 3’/^ chord 
point midway between the tips of the element. The strength of each bound 
vortex is chosen so that the total normal wash due to the oncoming free stream 
flow and to all horseshoe vortex elements is zero at all collocation points . 
Therefore, in the doublet lattice method the zero normal wash condition need 
be satisfied at only FT discrete points on the wing; nomaal wash may be present 
at all other points on the wing. 

The doublet lattice method does represent an efficient and versatile 
method of solving the inviscid flow about a wing of general planform and may 
be capable of giving accurate predictions of lift and moment coefficients. 
However, unless an extremely large number of rectangular sections and colloca- 
tion points are used, the procedure suffers from obvious deficiencies in its 
ability to determine the flow field in the immediate vicinity of the wing. 

The deficiency results from the zero normal wash boundary condition being 
satisfied only at isolated points (one per section) and, therefore, at all 
other points the correct physical boundary condition is not being satisfied. 

In the present inviscid calculation which' was used to obtain the inviscid 
pressure field required by the viscous procedure only a moderately dense grid 
was used and, consequently, the zero normalwash condition in the immediate 
vicinity of the wing tip was not well satisfied. A typical spanwise distribu- 
tion of the normalwash velocity is presented in Table I. The distribution 
was calculated by summing the contributions from each horseshoe vortex and 
the oncoming free stream flow. The distribution is taken at x/c = 0.325 and 
y/c = 0.02; the velocity is normalized by the oncoming free stream velocity. 
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and -O.Uc for the outboard and inboard transformations, respectively. It 
should be noted that the y-grid is ssnnmetric about the wing centerplane and 
the z-grid is symmetric about the wing tip. The locations of the grid points 
versus point number are presented in Pig. l6. 

In regard to the finite difference representation, in the present effort 
second derivatives are represented by three point central difference operators. 
First derivatives are represented by two different difference operators . In 
the spanwise direction a two point backwards difference operator is used. In 
the transverse y direction for points above the wing centerplane a two point 
backwards difference is used whereas for points below the wing centerplane a 
two point forwards difference is used. It should be noted that with this 
difference scheme , differencing is always done away from boundaries which in 
truth can only be set at infinity. 

The final item to the specified before proceeding to a discussion of 
results is the application of boundary conditions. During the course of the 
present effort several combinations of boundary conditions were investigated. 

One possible choice along all outer boundaries of the computational region is 
the specification of the inviscid function values. However, this specifica- 
tion should require the flow in the outer portions of the computational region 
to be free from both direct and indirect (viscous displacement) viscous effects. 
In the present computations, the regions in the vicinity of the upper and 
lower (y= constant) boundaries are not expected to be free from inner layer 
viscous displacement effects and, therefore, along these boundaries second 
derivatives of u, w, and P were set to zero. The most satisfactory boundary 
condition for the fourth variable, v, was obtained by setting the first- 
derivative, dv/3y, equal to the Inviscid value. On the left hand z-boundary 
(which cuts through the wing) the second derivatives of u, v and p are set 
to zero and the spanwise velocity, w is set equal to its inviscid value raodi- 
fied'by a boundary layer profile, f(y/S), in the vicinity of the wind surface. 
Finally, on the right hand z=constant boundary, which should be reasonably 
free from all viscous effects, P, u and v are set to their inviscid values 
and 9^/d^z is. set equal to zero. 

Once a calculation has begun to march downstream, the continuity equation 
and the three momentum equations are solved in finite difference form. How- 
ever, a special procedure is required to start the calculations. In this 
special starting procedure during the y-impliclt sweeps of the ADI procedure, 
the y-momentum equation is replaced by the condition that the pressure equal 
the inviscid .pressure . Likewise during the z-implicit sweeps the z-momentum 
equation is replaced by the inviscid pressure condition. This special starting 
procedure only is required on the first step of the calculation and serves to 
generate a viscous flow field compatible with the inviscid pressure 
distribution. 
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Ta'ble I - Calculated Inviscid Normalwash Velocity 
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Obviously, the zero normalwash condition is not being held in the immediate 
vicinity of the tip. Therefore, for the purpose of the viscous calculation, 
the tip was taken to be at z/c = -.1 and the z-location was then shifted by 
+0.1 so as to put the new assumed tip location back to z/c = 0. 

In addition to the problem of nonsatisfaction of the zero normalwash 
condition (which in principle could be overcome through a dense calculation 
grid), the doublet lattice method assumes the trailing vortices to lie in a 
plane which in the present inviscid calculation coincides with the plane of 
the wing. This restriction further constrains the flow since the vortices 
are not free to interact one with the other. Although such a constraint may 
not represent a significant problem in lift and moment coefficient calcula- 
tions, the constraint may lead to critical inaccuracies if the correct 
detailed flow field in the vicinity of the airfoil wing is the object of the 
investigation. Therefore, in summary the inviscid calculation procedure 
used has two drawbacks in so far as the wing tip flow field is concerned: 

(1) satisfaction of the zero’ normalwash condition only at specific points and 
(ii) constraint of the vortex trailing filaments to a specified plane. 

Predictions of the inviscid flow field obtained with the doublet lattice 
method are presented in Fig. 17-22. Figures 17-19 show contour plots of the 
streaiHwise, normal and spanwise (u, v, and w) velocity components, respectively, 
at x/c = 0.175- As can be seen in these figures the streamwise velocity above 
the wing is greater than that below the wing and the normal and spanwise flow 
velocities show flow around the wing. Flow around the wing can be deduced 
from the prediction of positive v at' the wing tip in conjunction with posi- 
tive w below the wing and 'negative w above the wing. In addition it should 
be noted that a strong positive v-velocity occurs just beyond the wing tip 
location for the inviscid flow calculation (z/c = 0.1). This reaches a maxi- 
mum value of approximately v = 0.2. Far from the wing the ratio of v/u 
approaches a value of approximately 0.1 which is consistent with a 0.1 radian 
incidence angle. It should be noted that the inviscid solution shows no evi- 
dence of a circular secondary flow pattern in the y-z plane. Above the wing 
the spanwise velocity w is consistently negative while below the wing it is 
consistently positive. The normal velocity is always positive. The lack of 
any circular inviscid flow pattern in the y-z plane is a result found at all 
stations examined and may be due to constraining the trailing vortices in a 
planar wake. Figures 20-22 show the three inviscid velocity components at the 
plane x/c = 0.375; the flow pattern is similar to that found in the upstream 
flow plane . 
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V iscous Solutiou - The inviscid flow field obtained from the doublet 
lattice procedure is used as input for the viscous calculation. The inviscid 
pressure field is imposed unchanged upon the streaniwise momentum equation and 
the inviscid velocity field is used to obtain both upstream initial conditions 
and some boundary conditions. Although the tip flow field in general will be 
turbulent, the present calculations were run as laminar at a chord Reynolds 
number of approximately 2000. It is expected that the tip flow field region 
will be qualitatively similar in the laminar and turbulent cases and, there- 
fore, a laminar prediction should serve to assess the three-dimensional 
viscous flow field generated by the code. Treatment of the turbulent problem 
would add the complication of hypothesizing a turbulence model and further- 
more would require a simultaneous assessment of the turbulence model and the 
basic three dimensional calculation. The separation of these two items in 
one calculation is difficult at best and, therefore, it is proper in this 
first assessment to avoid the complications of a turbulent flow. Turbulent 
flow can be considered at a later date by adding an appropriate turbulence 
model to the calculation procedure . 

The calculation was initiated at the station x/c = 0.175 by inputting 
the inviscid velocity field and adding a boundary layer correction in the 
vicinity of the wing. The computational grid used consisted of twenty-one 
points in the spanwise direction and forty-one points in the transverse 
direction. In the present calculation a boundary layer thickness, 6, was 
chosen to initiate the calculation and each velocity component, Uj, was 
scaled such that 


’Jj" INVISCID f(y/Sl (47) 

where y is the distance from the wing and f(y/S) is the Pohlhausen velocity 
profile. For the purpose of this preliminary assessment of the procedure, a 
convenient boundary layer thickness was simply assumed. Inboard of the wing 
tip (z < O) the boundary layer thickness was assumed to be c and 

outboard of the wing tip ( z > 0) Sg assumed to be .015 c. Although the 
initial plane flow field obtained by superiirg)osi-ng a Pohlhausen type profile 
upon an independently calculated inviscid flow does give a qualitatively 
reasonable set of initial plane conditions, the procedure may give profiles 
which change rapidly during the first few stations of the viscous flow cal- 
culation, This problem is likely to be particularly acute in the region of 
the airfoil tip. The generation of initial plane conditions which are in 
concert with both the inviscid pressure field and the viscous flow equations 
is a major problem which must be faced, in the three-dimensional viscous flow 
analysis . 
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The results of the viscous calculation are shown in Figs. 23-30* Jor 
the jr-z plane grid used, the calculation of one strearawise step required 
approximately 30 systena seconds of TMIVAC 1110 CAU time. The velocity field.s 
at the fifth streamwise station, x/c = 0.19, are presented in Figs. 23-25. 

The streamwise velocity, u, presented in Fig. 23 clearly shows the velocity 
boundary layers on both the upper and lower wing surfaces as well as at the 
wing tip. A detailed examination of the calculation shows the boundary 
layer thickness on the upper surface, 5]_/c, is approximately 0.08 and the 
boundary layer thickness at the tip is approximately 0.035. It should be 
noted that in order to define the boundary layer region, the scale in Fig. 

23 is considerably less than that in Fig. 17 . The major difference between 
the inviscid 'u' field of Fig, 17 and the' viscous 'u' field of Fig. 23 is 
the viscous requirement of a no slip boundary condition. Outside of the 
boundary layer region the inviscid and viscous predictions of the 'u’ com- 
ponent of velocity are very similar. The v and w components of velocity are^ 
shown in Figs . 2h and 25 . Again the scale of Figs . 2h and 25 are less than 
those of Fig. l8 and 19* However, a comparison of the detailed calculations 
of Figs. l8 and 2k show that the 'v* distribution is modified significantly 
by viscous effects. Both above and below the wing, V is modified by bound- 
ary layer displacement effects whereas in the vicinity of the wing tip, the 
rapid changes of 'v' with respect to z in the inviscid calculation are con- 
siderably softened by the viscous solution. Finally, the symmetry in the 
inviscid v.-component field above and below the wing is not present in the 
viscous case. This destruction of symmetry results from the boundary layer 
above the wing being subjected to an adverse pressure gradient Whereas 
that below the wing being, subjected to a favorable pressure gradient, -The 
major difference between viscous and inviscid solutions is found in the 
spanwise w-velocity calculations. In the inviscid calculation ’w’ is nega- 
tive at all points above the wing (Fig. 19), however, in the viscous cal- 
culation a thin layer of positive ’w’ appears immediately above the wing 
(pig. 25). At this stage since 'v' is positive above the wing, no circular 
flow pattern can be discerned. 

Viscous calculations at x/c = .35 presented in Figs. 26-30. The' 
region in the immediate vicinity of the tip is shown in detail in Figs . 

26-27. In particular, Fig. 26 shows the v-velocity component field and it 
should he noted that at this station negative values of 'v' have- appeared 
above the wing. Figure 27 shows that the region above the wing also contains 
both positive and negative values of 'w*. Therefore, a 'circular' type flow 
pattern has appeared with a center at the location where v = w = 0; this center 
point has the approximate coordinates y/c = .025, z/c = -.04, This 'circular' 
flow pattern does not contain the nearly circular streamlines expected in a 
tip vortex since w » v; furthermore, the vortex is relatively weak. However, 
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considering that the inviscid flow field prediction shows no evidence of a 
vortex C^’igs. 21 and 22 ) this emergence of a 'circular' flow pattern, due to 
viscous effects and pressure gradients generated by the viscous solution in 
the secondary flow plane, does indicate the potential capability of the 
three-dimensional viscous flow procedure in treating the airfoil tip problem. 

The overall flow field encompassing more than the -immediate tip region 
is shown in Figs. 28-30. When compared against the analagous calculations for 
x/c = .175 (Figs. 23 - 25 ), it is clear that the major difference is that at 
this latter station the boundary layers now encompass more of the flow field. 

In addition some spurious spanwise velocities have appeared at the outer 
calculation boundaries (ly/c|> . 3 ) which are beyond the locations shown in 
the figures. These velocities are unrealistic and indicate more effort must 
be devoted to the treatment of boundary conditions, however, they are confined 
to the outer part of the flow calculation region and, therefore, do not 
influence the immediate tip region shown in Figs. 26 and 27 . A table of 
calculated velocities in the vicinity of the wing tip at x/c = .354 is pre- 
sented in Appendix B. 


Summary and Conclusions for Viscous Procedure 

A three dimensional forward marching, viscous, subsonic flow calculation 
procedure has been applied to the airfoil tip problem. The procedure inte- 
grates a reduced set of Navler-Stokes equations in which (i) streamwise diffu- 
sion is neglected and (ii) the pressure gradient in the streamwise momentum 
equation is obtained from an external source. In the present effort the 
doublet lattice method was used to generate the streamwise pressure gradients 
and with the inviscid calculation grid used this method did not properly define 
the flow in the vicinity of the tip. Furthermore, the doublet lattice method 
constrains the trailing vortices to lie in a plane . 

These deficiencies in the inviscid pressure field make an assessment 
of the viscous procedure somewhat difficult. However, the viscous procedure 
has shown features expected in the flow such as the qualitatively proper 
development of the viscous boundary layers both on the upper and lower sur- 
faces of the airfoil as well as on the airfoil tip surface. The viscous cal- ' 
culation also shows the generally expected flow pattern around the tip from 
the lower to the upper surface of the airfoil. Finally, the viscous procedure 
does predict a 'circular' type flow pattern to appear above the airfoil suction 
surfacej this feature was completely lacking in the inviscid solution. There- 
fore, the procedure is promising in its application to the airfoil tip problem, 
however, the results presented must be regarded as preliminary. Further inves- 
tigations which would include improvements of the viscous as well as the-invis- 
cid procedure used to generate the required streamwise pressure gradients must 
be made before the present approach can be fully assessed. 
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APPENDIX A 


SOLIJTIOIT PROCEDURE 


Numerical Techniques 

As an outline of the particulars of the numerical method, the treatment 
of the continuity equation is considered, as this is the simplest equation, 
and yet this discussion will cover most aspects of the method. The flow region 
is discretized by grid points having equal spacings Ax, Ay, and Az. Provisions 
for nonuniform grid spacing will be introduced subsequently. The subscripts 
i, .j 'and superscript n are grid point indices associated with y, z, and x, 
respectively. Thus ^ ^ denotes 0 z.) where 0 can represent any of 

the dependent variablel. The subscripts are frequently omitted if clarity is 
preserved, so that 0 is equivalent of .. For convenience, the following 
shorthand difference operator notation il^used, for derivative difference 
formulas : 


2 

with analogous definitions for Here a parameter a has been intro- 

duced ( 0 ^ O' £ 1) so as to permit, continuous variation from backward to for- 
wardi differences. The standard central difference formula is recovered for 
a = 1/2. Throughout the following discussion, it is assumed that the solu- 
tion is known at and is desired at x 

Consider the continuity equation 


d(h^hjpu)y/ 3x + 3(hj hj/av)/ 3y + 3 (h|h2/>w)/ 3 2 = 0 (50) 

Equation (50) is differenced in the x or marching direction as follows: 


(h2h5/>u)"'^'-(h h />ur 


[3(h,hj^v)/ dy + d(h,h2/>w)/dz] 


Ax (l-t/3) 

+ (T^7 [3(h|h3/^/)/3y +3(h,h2^w)/a2] =0 


(50) 
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Here a parameter P has heen introduced so as to permit a variable centering of 
the scheme in the x direction. Eq.uation (51) produces a backward difference 
formulation for P = 0 and a Crauk-Nicolson formulation for p = 1. The ’depen- 
dent variables in Eq. (51) are linearized by expansion about the solution at 
x^ (Refs. 12 and 17). Here a first-order accurate linearization is used for 
the X derivative, and the result is 
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8,(h"| h^j/v) 

h"2/>''w'') *0 

(52) ■ 


After eliminating the pressure gradients in Eqs. b-d) via the equation of 
state, the procedure outlined above for the continuity equation can be employed 
to derive linear implicit different approximations analogous to Eq. (52) for 
the three momentum equations. The resulting difference approximations can be 
grouped by coordinate direction and written in the following compact linear 
matrix difference operator notation as 




Ax 


(53) 


where f is a column vector containing the dependent variables, p, u, v, w, and 
A is a square (4 x 4) matrix. D and are 4x4 matrices containing elements 
which are themselves spatial difference operators for the y and z directions, 
respectively. S is a column vector reserved for any source temns which may be 
present. The matrices A, Dy, and contain only quantities which are known 
from a computational viewpoint. Equation (53) is linear in 

The advantage in grouping the dependent variables by the direction of 
differencing is that n’jmerous ADI or splitting techniques are immediately 
available for reducing the multidimensional implicit equation (53) to a ' 
sequence of one-dimensional equations (e.g., Douglas & Gunn, Ref. 24; Yanenko, 
Ref. 21), and this permits efficient solution while retaining the favorable 
stability properties of the basic implicit scheme, the present application, 
however, the technique of splitting (Yanenko, Ref. 21) is being employed. 

Using the technique of splitting, Eq. (53) can be written as the following two- 
step scheme. 
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(54a) 


tf") tt Dy<I>*+S| 

AX ' ' 

_fL(<I)"+L4>*) e + S, (54b) 

Ax i z 

where §* is an intermediate resiilt having computational significance but no 
particular physical significance, and where -t- = S. 

As discussed in the main body of the report, for subsonic flow, an invis- 
cid solution is used as a first approximation to the pressure field. This is 
accomplished by setting p = p^ in the axial momentum equation, where is the 
inviscid pressure. In the axial momentxun equation, the pressure gradient is 
therefore replaced by 

dP - 

dx dx (55) 

The pressure gradient terms in the transverse mementum equations are not 
altered. 


Solutions of the Split Difference Equations 

The coupled set of linear Implicit difference equations arising along 
rows of grid points during each step of the ADI solution procedure, together 
with the prescribed boundary conditions, can be written in a form having the 
following matrix structure. 




For each grid point index i, 0^ is a column vector containing the dependent 
variables p, u, v, w. and axe square (V x 4 ) matrices containing 

the implicit difference coefficients, d^^ is a column vector containing only 
computationally known quantities. There are K + 1 grid points along the row 
under consideration. Difference approximations for the four governing equa- 
tions are associated with symbols having subscripts 1 through H- 1 , the sub- 
scripts 0 and W are associated with the boundary conditions , which may involve 
up to three grid points. Equation (56) represents 4 (iT + 1)- linear equation 
in 4 (n + 1.) dependent variables. Excluding the elements and Aj,, the 
matrix structure of Eq, (56) is block tridiagonal, and direct solution by 
standard block elimination techniques (cf., Isaacson 85 Keller, Eef. 22) is 
both straightforward and efficient. The precise scheme used here consisted of 
Gaussian elimination for a simple tridiagonal system (sometimes called the 
Thomas algorithm) but with elements of the tridiagonal matrix treated as 
square submatrices rather than as simple coefficients. The required inverses 
of diagonal suhmatrices were obtained by a Gauss-Jordan reduction. The 
additional operations necessary to include the nonblock-tridiagonal elements 
Cq and are easily incorporated provided the original block tridiagonal 
coding is carefully organized. 
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APPENDIX B - 


Calculated Velocities at X/C = 0.35^ 
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U - VEL 


\z/c= 

- 0.1042 - 

-0.0667 

-o.o 4 o 8 

-0.0234 

-0.0121 ■ 

. -o.oo 47 

y/c^^ 







O.iok-T 

1.0319 

1-0327 

1.0326 

1-0314 

1.0340 

i.o 426 

0.0833 

1.0131 

1.0155 

1.0160 

1.0121 

1.0091 

1.0257 

0.0656 

0.9610 

0.9654 

0.9668 

0.9603 

0.9560 

0.9940 

0.0510 

0.8451 

0.8509 

0.8533 

0.8480 

0.8588 

0.9327 

0.0390 

0,6844 

0.6913 

0.6952 

0.6949 

0.7263 

0.8266 

0.0293 

0.5150 

0,5223 

0.5272 

0.5325 

0.5768 

0.6849 

0.0214 

0.3632 

0.3697 

0.3747 

0.3832 

0,4291 

0.5310 

0.0151 

0.2382 

0.2431 

0.2473 

0.2562 

0.2960 

0.3844 

0.0100 

0.1384 

o.i 4 i 5 

0.1445 

0.1519 

0.1822 

0.2540 

0.0059 

0.0593 

0.0605 

0.0621 

0.0674 

0.0879 

0.1427 

0.0026 

0.0000 

0.0000 

0.0000 

0,0000 

0.0000 

o.oobo 

0.0000 

0.0000 

0,0000 

0.0000 

0.0000 

0.0000 

0.0000 

-0.0026 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-0.0059 

0.0736 

0,0736 

0.0740 

0.0764 

0.0889 

0.1295 

- 0.0100 ’ 

0.1656 

0.1657 

0.1663 

0.1682 

o.l 8 o 4 

0.2220 

-0.0151 

0.2737 

0.2742 

0,2750 

0,2764 

0.2874 

0.3274 

- 0.0214 

0.3981 

0.3991 

0.4003 

0 , 4014 . 

0.4102 

0.4462 

-0.0293 

0.5357 

0.5375 

0.5394 

o, 54 oi 

0.5465 

0.5773 

-0.0390 

0.6773 

0.6802 

0.6828 

0.6836 

0,6873 

0.7118 

-0.0510 

0 . 8 o 42 

. 0.8g81 

0.8116 

0.8128 

0.8146 

0.8316 

-0.0656 

0.8913 

0.8961 

0.9002 

0.9024 

o. 9 o 46 

0.9156 

-0.0833 

0.9280 

0.9331 

0.9375 

0.9410' 

0.9462 

0.9554 

- 0.1047 

0.9426 

0.9476 

0.9520 

0.9563 

0.9639 

0.9731 
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U - YEL (COTW’D) 


\z/c= 

0.0000 

0.0047 

y/c=\^ 



o.io 4 t 

1.0362 

1.0228 

0.0833 

1.0288 

1.0206 

0.0656 

1.0161 

1.0208 

0.0510 

0.9818 

1,0117 

0.0390 

0.9035 

0.9698 

0.0293 

0.7834 

0.8872 

0.0214 

o. 64 i 6 

0.7766 

0.0151 

.0.4990 

0.6591 

0.0100 

0.3670 

0.5501 

0.0059 

0.2496 

0.4551 

0.0026 

0.0006 

0.3363 

0.0000 

0.0000 

0.3210 

-0.0026 

0.0000 

0.3186 

-0.0059 

0.2190 

0.4002 

- 0.0100 

0.3066 

0.4594 

-0.0151 

. 0.4038 

0.5314 

- 0.0214 

0.5132 

0.6195 

-0.0293 

0.6341 

0.7206 

-0.0390 

0.7563 

0.8204 

-0.0510 

0.8615 

0.9000 

-0.0656 

0,9311 

0.9461 

-0.0833 

0.9611 

0.9606 

- 0.1047 

0.9735 

0.9657 


0.0120 

0.0234 

o.o 4 o 8 

1.0203 

1,0282 

1.0292 _ 

1.0227 

1.0320 

1.0341 

1.0255 

1,0372 

1..0420 

1 . 0213 

1 . 0398 ' 

1.0507 

0.9947 

1.0296 

1.0533 . 

0.9377 

0.9998 

1.0446 

0.8561 

0.9518 

1 . 024 i 

0.7654 

0.8942 

0.9954 

0,6797 

0,8369 

0.9644 

0.6056 

0.7862 

0.9356 

0,5447 

0,7503 

. 0,9131 

0.5305 

0.7381 

0.9033 

0,5263 

0,7321 

0.8972 

0.5559 

0.7421 

0.9003 

0.5991 

0.7675 

0..9094 

0.6532 

■ 0.7988 

0.9195 

0.7209 

0.8385 

0.9332 

0.7992 

0.8851 

0.9505 

0.8734 

0.9268 

0.9644 

0.9281 

0.9542 

0.9713 

0.9555 

0.9654 

0.972I 

0.9606 

0.9656 

0.9701 

0.9602 

.0.9642 

0,9686 
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V - VEL 



-0.1042 

-0.0667 

-0.0408 

-0.0234 

-0.0121 

-0.004? 

Y/cX 







0.104? 

■ 0,01780 

0.01341 

0.01164 

0.02088 

o.o4?6o 

0.06144 

0.0833 

0.01357 

0.00988 

0.00859 

0.01766 

0.04413 

0.06093 

0.0656 

0.00828 

0.00535 

0.00451 

0.01303 

0.03814 

0.05716 

0.0510 

'0.00338 

6.00118 . 

0.00073 

0.00838 

0.03065 

0.04983 

0.0390 

-0.00017 

'-0.00176 

-0,00191 

0.00463 

0.02304 

0.04033 

0.0293 

-0.00200 

-0.00313 

-0.00312 

0.00211 

0.01621 

0.03022 

0.0214 

-0.00241 . 

-0.00318 

-0.00309 

0.00070 

0.01053 

0.02045 

0,0151 

-0.00198 

-0 .'00245 

-0.00235 

0.00001 

0.00580 

0.01124 

0.0100 

-0.00115 

-0.00138 

-0.00133 

-0.00034 

0.00176 

0.00279 

0.0059 

-0.00015 

-0.00019 

-0.00024 

-0.00040 

-0.00li2 

-0.00350 

0.0026 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.0000 ■ 

0.00000 

0.00000 

0., 00000 

0.00000 

‘o'. 00000 

0.00000 

-G.0026 

0.00000 

0.00000 

0.00000 

0 .. 00000 

0.00000 

0,00000 

-0.0059 

-0.00024 

-0.00022 

-0.00018 

0.00006 

0.00125 

0.00510 

-0.0100 

0.00010 

0.00020 

0.00015 

-0.00040 

-0.00152 

-0; 00070 

-0.0151 

-0.00004 

0.00019 

0.00014 

-0,00106 

-0.004?2 

-0.00773 

-0.0214 

-0.00082 

-0.00042 

-0.00042 

-0.00208 

-0.00784 

-0.01435 

-0.0293 

-0.00239 

-0.00181 

-0., 00174 

^0.00376 

-0.01119 

-0.02044 

-0.0390 

-0.00488 

-o.oo 4 i 3 

-o,oo 4 o 3 

-0.00639 

-0.01516' 

-0.02631 

-0.0510 

-0.00841 

-0.00753 ■ 

-0.00742 

-0.01013 

-0.01990 

-0.03189 

-0.0656 

-0.01275 

-0.01180 

-0.01172 

-0.01478 

-0.02510 

-0.03655 

-0.0833 

-0.01704 

-0.01610 

-0.01611 

-0.01947 

-0.02978 

-0.03918 

-0.104? 

-0.02000 

-0.01918 

-0.01935 

-0.02296 

-0.03263 

-0.03883 
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V - VEL (COKTINUED) 


^\Z/C = 
Y/C 

0.0000 

0.004? 

0.0121 

0.0234 

O.o4o8 

0.104? 

0.03993 

-0.01614 

0.01755 

0.04487 

0. -05134 

' 0.0833 

0.03995 

-0.02577 

0.01439 

0.04417 

. 0.05171 

0.0656 

0.03864 

-0.03358 

0.0104? 

0.04273 

0.05183 

0.0510 

0.03495 

-0.03880 

0.00544 

0.03999 

0.. 05120 

0.0390 

0.02852 

-0.04l4l 

0.00016 

0.03617 

0.04969 

0.0293 

0.02049 

-0.04126 

-0.00424 

0:03206 

■ 0.04771 

0.0214 

0.01206 

-0.03857 

-0.00673 

0.02873 

0,04594 

0.0151 

0.00389 

-0.03395 

rO. 00669 

0.02704 

0.o4504 

0.0100 

-0.00362 

-O'. 02806 

-0.00402 

,0.02734 

0.04528 

0.0059 

-0.00898 

-0.02081 

0.00129 

0.02966 

0.04664 

0.0026 

0.00000 

-0.00244 

0.01307 

0.03429 

0.64821 

0.0000 

0.00000 

0. 00776 

0.02039 

0.03772 

0.04964 

-0.0026 

0.00000 

0.01771 

0.02773 

0.04123 

0.05113 

-0.0059' 

0.01403 

0.03363 

0.03836 

0.04613. 

0.05328 

-0.0100 

0.00881 

0.03887 

0.04590 

0.05146 

0.05627 

-0.0151 

0.00211 

0.04292' 

0.05205 

0.05595 

0.05894 

-0.0214 

-o.oo4i4 

0-. 04?04 

0.05762 

0.06011 

0.06152' 

-0.0293 

-0.00950 

0., 05148 

0.06282 

0.06401 

0.06403 

-0.0390 

-0.01410 

0.05537- 

■ 0.. 06706 

0.06720 

0.06614 

-0.0510 

-0.017T1 

0.-05773 

0.06993 

0.06955 

0.06780 

-0.0656 

-0.01966 

0.05798 

0.07122 - 

0.07106 

0.06908- 

-0.0833 

-0.01925 

0.05610 

0.07091 

0.07174 

0.06999 

-0.104? 

-0.01641 

0.05265 

0.06922 

0.07166 

0.07052' 


i ■■ n < tl xKj’ I- ' - 

' ‘ u*I :v« 
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W - TEL 


\z/c = 

- 0.1042 

-0.0667 

- 0.0408 

-0.0234 

- 0.0121 

-0.0047 

0.1047 

-0,0826 

‘ -0.0927 

- 0.1000 

- 0.1046 

-0.1072 

-0.1086 

0.0833 

-0.0871 

-0.0974 

- 0.1048 

-0.1093 

-0.1115 

-0.1130 

0.0656 

-0.0815 

-0.0901 

-0.0962 

-0.0999 

-0.1015 

-0.1037 

0.0510 

-0.0633 

-0.0685 

-0.0722 

-0.0745 

-0.0759 

-0.0803 

0.0390 

- 0.0406 

- 0.0421 

- 0.0435 

- 0.0444 

- 0.0461 

-0.0526 

0.0293 

-0.0180 

-0.0163 

-0.0155 

-0.0151 

-0.0172 

- 0.0246 

0.0214 

0.0030 

0.0075 

0.0104 

0.0120 

0.0100 

0.0024 

0.0151 

0.0233 

0.0305 

0.0356 

0.0383 

- 0.0367 

0.0294 

0.0100 

0.0469 

0.0574 

0.0646 

0.0687 

0.0678 

0.0605 

0.0059 

0.0803 

0.0953 

0.1056 

0.1113 

0.1112 

0.1034 

0.0026 

0.0000 

0.0000 

0.0000 

0.0000 

o.oobo 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0-0000 

0.0000 

-0.0026 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0,0000 

-0.0059 

0.0963 

0.1119 

0.1220 

0.1278 

0.1288 

0.1224 

- 0.0100 

0.0748 

0.0864 

0.0942 

0.0988 

0.1001 

0.0959 

-0.0151 

0.0671 

0.0762 

0.0824 

0.0861 

0.0873 

0.0846 

- 0.0214 

0.0657 

0.0730 

0.0779 

0.0808 

0.0818 

0.0800 

-0.0293 

0.0639 

0.0695 

0.0730 

0.0751 

0.0759 

0.0745 

-0.0390 

0.0590 

0.0623 

0.0644 

0.0655 

0.0658 

0.0646 

-0.0510 

0.0507 

0.0517 

0.0521 

0.0521 

0.0520 

0.0510 

-0.0656 

0.0399 

0.0387 

0.0375 

0.0366 

0.0360 

0.0351 

-0.0833 

0.0284 

0.0254 

0.0230 

0.0213 

0.0201 

0.0192 

- 0.1047 

0.0184 

0.0141 

0.0107 

0.0084 

0.0069 

0.0058 
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W - VBL (CONT'D) 


= 

0.0000 

0.0 o 47 

0.0121 

0.0234 

o.o4o8 

o.io 47 

-0.1093 

-0.1102 

-0.1093 

-0,1088 

-0.1086 

0.0833 - 

-0.1144 

-0.1157 

-0.1158 

-0.1163 

-0.1174 

0-0656 

-0,1068 

-0.1096 

-0.1111 

-0.1136 

-0.1170 

0.0510 

-0.0857 

-0.0909 

-0.0937 

-0.0932 

-0.1047 

0.0390 

-0.0597 

-0.0670 

-0.0698 

-0.0755 

-0,0835 

0.0293 

-0.0327 

-o.o 4 i 7 

-0.0430 

-0.0483 

-0.0566 

0.02 i 4 

-0.0063 

-0.0171 

-0.0159 

-0.0195 

-0.0266 

0.0151 

0.0195 

0.0057 

0,0107 

0.0098 

0.0051 

0.0100 

0.0484 

0.0283 

0.0389 

0.0423 

o.o 4 o 9 

0.0059 

0.0866 

0.0543 

0.0740 

0.0838 

0.0860 

0.0026 

0.0000 

-0.0001 

0.0531 

0.0800 

0.0923 

0.0000 

0.0000 

0.0002 

0.0519 

0.0792 

0.0930 

-0.0026 

0.0000 

0.0005 

0.0527 

0.0809 

0.0959 

-0.0059 

0.1063 

0.0750 

0.0962 

0.1109 

0.1190 

-0.0100 

0.0857 

0.0676 

0.0813 

0,0920 

0.0995 

-0.0151 

0.0780 

0.0676 

0.0774 

0.0860 

0.0927 

-0.0214 

0,0756 

0.0692 

0.0766 

0.0835 

0.0890 

-0.0293 

0.0713 

0.0670 

0,0725 

0.0774 

0,0808 

-0.0390 

0.0622 

0.0592 

0.0628 

0.0657 

0.0668 

-0.0510 

0.0492 

0.0471 

0.0491 

0.0503 

o,o 499 

-0.0656 

0.0339 

0.0326 

0.0336 

0.0338 

0-:0329 

-0.0833 

0.0185 

0.0176 

0.0183 

0.0181 

0.0174 

-0.1047 

0.0053 

0.0044 

0.0052 

0.0050 

0.0044 
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Figure 1. — Schematic of vortex formation at tip of wing with finite thickness. 
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Figure 2. — Schematic of vortex formation at tip of thin flat plate 
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Figure 4. — Geometry of encounter of infinite wing with gust wave fronts 
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a) RELATIVE MACH NUMBER Wl 




Figure 5. — Geometric definitions 
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Figure 6. Comparison of lift transfer functions 
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Figure 8. — Real and imaginaty part of pressure difference coefficient for IVI = 0.3 


60 




REAL AND IMAGINARY PART OF AC 



X/C 


a 


O 


b) VARIATION WITH REDUCED FREQUENCY 




Figure 9, — Real and imaginary part of-pressiire difference coefficient for M = 0.6 
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Figure 10. Effect ofspanwiseppsition-on steady chord wise, pressure distribution 
for an.unswept semi-infinite wing 
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Figure 13. — Unsteady response of flat plate semi-infinite wing at zero angle of attack 
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Figure 14. — Pressure distribution in vicinity of tip of semi-infinite wing 
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Figure 15. — Viscous flow coordinate system. 
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Figure 16. — Grid point locations. 
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Figure 18, — Inviscid flow calculation 
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Figure 19. — Inviscid flow calculation. 
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Figure 20. Inviscid flow calculation, 
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Figure 21 
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Figure 23. — Viscous flow calculation. 
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Figure 25. — Viscous fiow calculation. 



SPANWISE DISTANCE, Z/C 


Figure 26. — Viscous flow calculation, 
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Figure TCI. — Viscous flow calculation. 
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Figure 29. — Viscous flow calculation 
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